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Abstract 

A computational procedure that allows the detection of a new type of high- 
dimensional chaotic saddle in Hamiltonian systems with three degrees of 
freedom is presented. The chaotic saddle is associated with a so-called nor- 
mally hyperbolic invariant manifold (NHIM). The procedure allows to com- 
pute appropriate homoclinic orbits to the NHIM from which we can infer 
the existence a chaotic saddle. NHIMs control the phase space transport 
across an equilibrium point of saddle-centre-... -centre stability type, which 
is a fundamental mechanism for chemical reactions, capture and escape, 
scattering, and, more generally, "transformation" in many different areas 
of physics. Consequently, the presented methods and results are of broad 
interest. The procedure is illustrated for the spatial Hill's problem which is 
a well known model in celestial mechanics and which gained much interest 
e.g. in the study of the formation of binaries in the Kuiper belt. 

PACS: 05.45.Jn, 45.50.Pk 
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Introduction 



Chaotic saddles are the saddle-type invariant Cantor sets associated with a 
horseshoe construction ^1^. They play a central role in many complex dy- 
namical phenomena, e.g., the existence of supertransients j2] and the fractal 
structm'e of chaotic scattering 0]. The chaotic saddles constructed to date 
have been related to either homoclinic orbits associated with hyperbolic pe- 
riodic orbits or special types of equilibria. In this letter we are concerned 
with a fundamentally new type of high dimensional chaotic saddle structure 
that can occur only in Hamiltonian systems with three or more degrees of 
freedom (DOF), and we describe a computational method for detecting it 
which also illustrates its geometrical complexity. 

The chaotic saddle is associated with a normally hyperbolic invariant 
manifold (NHIM) [3]. The physical significance of NHIMs arises from the 
fact that their stable and unstable manifolds control the phase space trans- 
port across equilibria of saddle-centre-... -centre stability type. This is is not 
only the fundamental mechanism for the evolution from reactants to prod- 
ucts in chemical reaction, but also for "transformations" in general in a large, 
and diverse, number of applications as e.g. ionisation problems in atomic 
physics [3j , rearrangements of clusters |6| , cosmology [7| , and solid state and 
semi-conductor physics jHlini- Though it had been recognised that it is im- 
portant to understand the dynamics near saddle-centre-... -centre equilibria 
it was only recently that new developments in dynamical systems theory 
offered the theoretical framework and computing power offered the means 
to study the phase space structure near saddle-centre-... -centres for systems 
with 3 or more DOF [TOl HH dUS] . Besides the NHIM and its stable and 
unstable manifolds it is now possible to compute a codimension 1 submani- 
fold of the energy surface, the so-called transition state, which is transverse 
to the Hamiltonian flow and locally divides the energy surface into two dis- 
joint components. Locally, the transition state provides the only means of 
passing from one phase space region (associated with "reactants" ) to another 
phase space region (associated with "products"), i.e. trajectories must cross 
the transition state in order to "react" . The transversality of the transition 
state to the Hamiltonian flow is essential for rate calculations as it solves 
the problem of locally recrossing trajectories. 

In this letter we go beyond this local result and consider more global is- 
sues associated with the dynamics related to the transition state and NHIM. 
In particular, we describe a computational method for determining the ex- 
istence of homoclinic and heteroclinic connections to NHIMs. Using recent 
results in ^3], the existence of certain types of these homoclinic orbits allows 
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us to infer the existence of a new type of chaotic saddle (in fact, a Cantor set 
of chaotic saddles). Because of the ubiquity of saddle-centre-... -centre type 
equilibria in applications (as described above), we expect that the methods 
and results presented here, which are applicable to three and more DOF, to 
be of broad interest. 

The physical system we choose to illustrate our method is one from celes- 
tial mechanics; the 3-d Hill's equations [E]. Advances in detector technology 
have opened up new frontiers in celestial mechanics with the discovery of 
trans-Neptunian objects and binary systems in the Kuiper belt. The calcu- 
lation of capture probabilities requires the study of the 3-d Hill's equations 
(rather than the thoroughly studied 2-d case) because many of the capture 
events occur from high inclination (see e.g. the recent work by Goldreich, 
Lithwick and Sari [T6]). 

Hill's Problem and the Phase Space Structure Near 
Saddle-centre-centre Equilibria 

The circular restricted three body problem (CRTBP) models the motion of 
a tiny particle under the gravitational influence of one (large) primary mass 
and one (smaller) secondary mass both in circular orbits about their common 
centre of mass |17| . Hill's problem is a limit version of the CRTBP which 
describes the motion of the particle in a neighbourhood of the secondary 
mass. Hill's equations can be derived from the following Hamiltonian in 
dimensionless units, 



two equilibria, which are denoted traditionally by Li and L2, and that the 
matrices associated with linearising Hamilton's equations about each equi- 
librium have a pair of real eigenvalues and two pairs of pure imaginary 
complex conjugate eigenvalues. This means that Li and L2 are equilibria of 
saddle-centre-centre type. 

A detailed theory for phase space transport near saddle-centre-centre 
equilibrium points has been developed in recent years ^1 ^1 • For 
energies slightly above that of the saddle-centre-centre equilibrium point, 
on the 5-dimensional energy surface there exists an invariant 3-dimensional 
sphere of saddle stability type. This 3-sphere is significant for two rea- 
sons: 



H = 





It is well known that Hill's equations have 
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• It is the "equator" of a 4-dimensional sphere cahed the transition state. 
Except for the equator (which is an invariant manifold), the transition 
state is locally a "surface of no return" in the sense that all trajectories 
that start in a neighbourhood of the transition state in the energy 
surface must cross the transition state and exit the neighbourhood (in 
the appropriate direction in time). For energies "sufficiently close" to 
the energy of the saddle-centre-centre equilibrium point, the transition 
state satisfies the bottleneck property. This means that, the energy 
surface locally has the geometrical structure of S"^ x / (i.e., 4-sphere 
X interval) and the transition state divides the energy surface into 
two disjoint components. Moreover, the only way a trajectory can 
pass from one component of the energy surface to the other is to pass 
through the transition state. 

• The 3-sphere is a normally hyperbolic invariant manifold j4j (NHIM), 
which means that the expansion and contraction rates of the dynamics 
on the 3-sphere dominate those transverse to it. Just like a "saddle 
point" it therefore has stable and unstable manifolds. In this case the 
stable and unstable manifolds are 4-dimensional, having the structure 
of spherical cylinders, x M. Hence, they are of one less dimension 
than the energy surface and act as "separatrices" ; they "enclose" a 
volume of the energy surface. Their key dynamical significance is that 
the only way that trajectories can pass through the transition state 
is if they are inside the region of the energy surface enclosed by the 
stable and unstable spherical cylinders. 

These phase space structures can be realised through a procedure based 
on Poincare-Birkhoff normalisation by which explicit formulae for the NHIM, 
its stable an unstable manifolds, and the transition state are given jl3j in 
"normal form coordinates". The phase space structures are then mapped 
back into the original coordinates by the inverse of the normal form (NF) 
transformation. 

Computation and Visualisation of the Phase Space 
Structures near Li and L2 

We will visualise the phase space structures near Hill's equilibria Li and 
L2 as their projections to configuration space. In understanding the result, 
and its implications, it is useful to recall that the level sets of the effective 
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potential energy 



y = H - ^\{px + y,Py - x,pz)\^ = -^x^ + - ^, 

the so-called zero velocity surfaces, confine the motion in configuration space. 

In figure Wi&) we show the zero velocity surface (ZVS) V = E for an 
energy E = —4.4 which is 0.1 above the energy of the equilibrium points. 
The ZVS encloses the energetically-allowed volume of configuration space. 
For Hill's equations its shape has two "bottlenecks" associated with Li and 
L2 which divide the allowed volume into three regions: a bounded region 
about the origin (for which — 1 < x < +1 (A)) and two unbounded regions 
(for which x < —1 (B) and x > +1 (C), respectively). The phase space 
structures near Li and L2 regulate the transport through the bottlenecks 
between these regions. 

The NFs about Li and L2 are related by symmetry. It is thus sufficient 
to compute explicitly only the NF about Li. As a result of the two complex 
eigenvalues associated with the equilibria being rationally independent, the 
NF, to any desired finite order of computation, is completely integrable with 
integrals given by T = piqi, Ji = \ [qf + pf) , i = 2,3. The NF Hamilto- 
nian can be written solely as a function of the integrals, H = H(Z, J2, J3), 
and Hamilton's equations decouple into the product of independent linear 
systems 

OH dH 
iQi,Pi) = -Qjili^ -Pi)^ {Qi,Pi) = -QjiPi^ ^ = 2,3. 

The transition state 4-sphere is given by qi = pi- From the equations of 
motion it is easy to see that it is a "surface of no return". Its equator 
Pi = qi = is the NHIM (3-sphere). It has 4-dimensional stable {qi = 0) 
and unstable (pi = 0) manifolds. In the energy surface volume enclosed 
by the stable and unstable manifolds X is positive; outside I is negative. 
The NHIM has a special structure. It is foliated by a one parameter family 
of invariant 2-tori (the "Hopf fibration") jlUL 1111 113j . These tori can be 
parametrised e.g. by the integral J2 where J3 is then given implicitly by 
energy conservation and 1 = 0. At the minimal and maximal values of J2 
the 2-tori degenerate to periodic orbits, the so-called "Lyapunov orbits". 
The 2-tori and the periodic orbits have 3-dimensional and 2-dimensional 
stable and unstable manifolds, respectively, which are contained in the 4- 
dimensional stable and unstable manifolds of the NHIM. These geometrical 
considerations will play an important role described below. 



Fig. 1 
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We perform the NF computation to order 18 using the computer alge- 
bra system Mathematica. As a result the NF Hamiltonian is a sum over 
219 multivariate monomials in (I, J2, J3). Each component of the mapping 
between NF coordinates and the original phase space coordinates involves 
sums over about 50 000 multivariate monomials. The resulting transition 
state and NHIM are shown in figure Q^b). Notice that the transition state 
completely blocks the "bottleneck" in the ZVS in configuration space; it also 
blocks the bottleneck in the level set of the Hamiltonian in phase space. 

Computation of Homoclinic and Heteroclinic Con- 
nections between the NHIMs near Li and L2 

The "saddle integral" , I = piqi, plays the key role in our numerical approach 
to detecting orbits connecting the same NHIM, or orbits connecting different 
NHIMs. Taking into account the Hopf fibration of the NHIM our procedure 
consists of the following four steps: 

• In the NF coordinates, choose an invariant 2-torus on the NHIM and 
seed a mesh covering the torus with initial conditions. Displace these 
initial conditions slightly in the direction of the unstable manifold of 
this torus {pi = 0, qi = e). 

• Map the initial conditions back into the original coordinates using the 
NF transformation. 

• Integrate the initial conditions forward in time using Hill's equations. 
Since they are in the unstable manifold they will leave the neigh- 
bourhood in which the NF transformation is valid (which is why we 
integrated them in the original coordinates). 

• Check if a trajectory returns to the neighbourhood of Li or L2 where 
the NF is valid. If so, map it into the NF coordinates and check the 
value of its saddle integral. If the saddle integral is zero, the trajec- 
tory must be on the stable manifold of the respective NHIM (since 
the trajectory is already on the unstable manifold and two unstable 
manifolds cannot intersect). 

The procedure can be understood as a shooting method between the -^^S- 
X-fibres of the locally valid NFs. We apply the method to Hill's problem by 
displacing (with e = 10"'*) meshes of initial conditions on invariant 2-tori of 
the NHIM near Li along the respective unstable manifold branches which 
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are directed towards the bounded region about the origin, see figure^ These 
initial conditions are propagated by integrating Hill's equations where the 
singularity at the origin is taken care of by the Kustaanheimo-Stiefel regu- 
larisation. If the resulting trajectory reenters or enters the neighbourhood 
of validity of the NF about Li or L2, respectively, we check the value of the 
integral X. A positive X means that the trajectory will cross the respec- 
tive transition state which leads to an exit to the outside of the bounded 
region about the origin. In this case the integration is stopped. A negative 
value of X means that the trajectory does not exit on this approach of Li 
or L2 and the integration is continued until the trajectory reaches a validity 
neighbourhood of the NF with positive X. It is to be noted that along the 
part of a trajectory which traverses the validity neighbourhood of the NF 
the values of the integrals (I, J2, J3) are conserved to 12 digits and more, 
demonstrating the high accuracy of the NF. 

We illustrate the results in figure |2l where contours of the integral X\^^ 
(light blue/dark blue; exit through the transition state near Li) and X\^^ 
(yellow/red; exit through the transition state near L2) are shown plotted on 
the corresponding torus of the NHIM. At the boundary of a light blue/dark 
blue region goes to zero and points on the boundary correspond to 
homoclinic orbits which connect the NHIM near Li to itself. Similarly, 
boundaries of yellow/red regions represent heteroclinic connections between 
the NHIMs near Li and L2. 

Though the individual regions in figure themselves are regular in the 
sense that they each have a smooth boundary, the disposition of the regions 
is very intricate. We illustrate this for initial conditions on the 2-dimensional 
unstable manifold of the Lyapunov periodic orbit near Li which has J3 = 
and in configuration space lies in the (x, y)-plane, see figure 13 Between 
each two arcs which extend from two consecutive zeros of Z as a function 
of the angle conjugate to J2 there is an infinity of further arcs. Each zero 
represents an orbit which in time is backward asymptotic to the Lyapunov 
periodic orbit and forward asymptotic to the NHIM near Li (blue) or the 
NHIM near L2 (red). The result is a self-similar structure well known from 
classical scattering systems. 

Homoclinic Connections and Chaotic Saddles 

While orbits homoclinic to normally hyperbolic invariant tori have been 
studied (Q^), tori cannot be normally hyperbolic in Hamiltonian systems 
jl9j . The NHIM is normally hyperbolic, yet there are no theorems describ- 
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ing the dynamics associated with orbits homochnic to a normally hyperbolic 
invariant sphere (part of the difficulty here comes from the fact that a sphere 
cannot be described by a single coordinate chart, but see 'T._ for numerical 
evidence that this could be an important mechanism for chaos). An im- 
portant piece of this problem has recently been solved by Cresson jl4j who 
proved that if the stable and unstable manifolds of a torus in the NHIM 
intersect transversely, then there exists a (uniformly) hyperbolic invariant 
Cantor set on which the dynamics is conjugate to a shift map, i.e., a chaotic 
saddle. The stable and unstable manifolds of a 2-torus are 3-dimensional 
in the 5-dimensional energy surface. A transverse intersection is necessarily 
1-dimensional, i.e., a trajectory. 

When the effect of the neglected terms in the normal form expansion 
(i.e., the non-normalised "tail" of the expansion) are included we expect, by 
KAM theory, that a Cantor set of non-resonant tori in the Hopf fibration will 
persist. Homoclinic intersections of the stable and unstable manifolds of a 
given 2-torus (rather than more general connections between the manifolds 
of such a torus and those of the entire NHIM) can be detected easily by 
our numerical procedure described above. We therefore not only check the 
value the integral X upon entering the neighbourhood of validity of the same 
NF that we started from, but also the difference between the centre integral 
J2 of the 2-torus we started from and the centre integral J2 in the entered 
validity neighbourhood. If both J2 — J2 and I vanish we have a homoclinic 
orbit to the 2-torus under consideration. The result of such a procedure is 
shown in figure|Ha) where homoclinic connections appear as the intersection 
points of the zero contours of J2 — J2 (sharp boundaries between dark blue 
and light blue within the blue regions) and those of T (boundaries of the 
blue regions themselves). Figure Efb) shows the example of a homoclinic 
orbit which corresponds to one of the prominent intersection points marked 
by a white dot in the figure. In fact there are many more intersection points 
located in the finer structures of figure |lf a). There is a tendency that the 
homoclinic orbits become more complicated in the finer regions. 



Conclusions 

In this letter we have described a numerical method for detecting a Cantor 
set of chaotic saddles in Hamiltonian systems with three DOF. We have 
illustrated it for the spatial Hill's problem. The implications for this prob- 
lem should be of current interest in celestial mechanics. Recently it has 
been shown that chaos plays an important role in the energetics of cap- 
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ture [20], and this Cantor set of chaotic saddles should be central to this 
process. Chaotic scattering in systems with three or more DOF is poorly 
understood, and our results and methods should provide a window into this 
subject. Towards this end, our methods detect the chaotic saddles by de- 
tecting appropriate homoclinic orbits. 

Concerning transition state theory, as described by Truhlar j^, there 
are two types of trajectories that exhibit multiple recrossings of the transi- 
tion state: local and global. As described in the introduction, the correct 
choice of transition state solves the local recrossing problem. However, char- 
acterising the global recrossing problem, i.e. the question of whether there 
are trajectories, which, after crossing the transition state and leaving its 
neighbourhood, return to the transition state and cross it again, is more 
difficult. Under additional assumptions on the geometry of the homoclinic 
or heteroclinic orbits, the chaotic saddles consist of trajectories that enter 
and leave a neigborhood of the transition state, recrossing it infinitely often. 
This will be the subject of a future publication. 

It would also be of some interest to actually compute and visualise the 
chaotic saddles themselves. A recently developed numerical method j21j 
may play a role for this purpose. Finally, we remark that our results are 
not limited to 3 DOF. The general theory developed in ^3 ^3 El Ej 
applies to Hamiltonian systems with arbitrary DOF. The numerical method 
generalises easily since the saddle integral is a scalar function regardless of 
the number of DOF. 
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Figure Captions 



Figure^ (a) Cut-away of the zero velocity surface (blue/green), (b) Tran- 
sition state at L2 (red) and NHIM at Li (white). The energy is -4.4. 

Figure [2J The top panels highlight individual 2-tori in the Hopf fibration 
of the NHIM near Li and shows contours of values of the saddle integral 
and on the 2-tori. The 2-tori are for J2 = nJ2max/4, n = 1,2,3, where 
•^2 max is the maximum J2 on the NHIM. n = and n = 4 correspond to the 
two Lyapunov periodic orbits. The 2-tori are parametrised by the angles 
and 03 conjugate to J2 and J3. {E = —4.4.) For clarity the bottom panels 
show the tori in the covering space. 

Figure El Homoclinic and heteroclinic orbits between the planar Lyapunov 
orbit near Li and the NHIMs near Li (zeros of blue arcs) and L2 (zeros of 
red arcs). 

Figure 12 (a) Points of homoclinic connections on the 2-torus J2 = J2max/2 
(see figure 121) • (b) A sample homoclinic trajectory in configuration space. 
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